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ABSTRACT 
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A change-of-basi s transformation  between  the  trigonometric  Fourier 
Series  and  the  Legendre  polynomial  series  is  derived.  This  transformation 
leads  to  two  immediate  results:  one  is  a fast  algorithm  for  computing  the 

Legendre  polynomial  coefficients  for  a function;  the  second  is  a numerically 
stable  computation  algorithm  for  evaluating  jf,(z),  the  n^*^  order  spherical 
Bessel  function.  In  addition,  these  results  lead  to  a simple  recursive 
digital  network  to  be  used  in  the  implementation  of  adaptive  nonlinear 
functions. 

I.  INTRODUCTION 

Two  common  ways  of  representing  functions  have  been  polynomial  expan- 
sions and  trigonometric  expansions.  In  much  of  electrical  engineering  the 
trigonometric  expansion  has  useful  Interpretations  and  has  dominated  over 
the  generalized  Fourier  series  expansions  in  applications.  However,  many 
functions  are  readily  expressed  in  terms  of  polynomials.  For  example,  in  a 
recent  paper,  we  have  presented  an  interesting  representation  for  bandl imi ted 
functions  using  the  Legendre  polynomial  expansion  [1].  In  this  paper  we 
derive  a simple  linear  transformation  which  maps  the  polynomial  representa- 
tion into  a trigonometric  representation.  Also,  we  derive  the  inverse 
transformation  which  maps  a trigonometric  expansion  to  a polynomial  expan- 
sion. 

The  inverse  transformation  has  enabled  us  to  develop  a fast  algorithm 
for  the  computation  of  the  Legendre  polynomial  coefficients  for' any  Ln 
[-ir,ir]  function.  The  algorithm  utilizes  the  Fast  Fourier  Transform  (FFT) 
to  compute  the  Fourier  series  coefficients  and  then  multiplies  the  vector 
of  coefficients  by  a linear  matrix  transformation  to  compute  the  vector  of 
polynomial  coefficients.  This  approach  can  offer  a considerable  saving  in 
computation  time  over  the  standard  integral  formula  for  computing  these 
polynomial  coefficients. 

Section  II  contains  the  derivation  for  the  elements  of  the  transforma- 
tion matrix  and  its  inverse.  In  section  1 1 1 we  discuss  how  these  results 
lead  to  the  design  of  a simple  recursive  digital  filter  that  may  be  employed 
to  lnq>iement  almost  any  zero  memory  non-linear  operation  (ZNL).  In  section 
IV  we  discuss  several  computational  considerations,  and  in  section  V we 
present  an  example. 

1 1 . DEVELOPMENT 

Assume  H(x)  is  a polynomial  defined  on  I-ir,ir].  By  writing  H(x)  In 
Legendre  polynomials,  we  can  use  results  developed  in  an  earlier  paper  to 
easily  write  the  Fourier  series  expansion  for  H(x)  that  converges  In  t-w,irl. 

PfLe^^nJt^d  at  thz  PouAtzznth  Annual  AtleAton  Con^eAcncz  on  CviciUt  and 
Sf/Atem  Thzofiy,  1976;  to  be  pubtUhe.d  in  thz  P/toceedcngA  the  Conitfience. 


Assume  that 
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"f-'  ■ ' “n  %'ri- 

n-0 


n < X < n 


where  Legendre  polynomial.  To  determine  the  Fourier  series 

representation  for  H(x),  it  is  sufficient  to  know  the  Fourier  series  expan- 
sion for  P„(— ).  From  [1]  we  know  that 

**  n 

P„(7)  - r KD"  j„(nm)]  e’*"”‘  . (2) 

m“-« 

By  Eqs.  (1)  and  (2),  we  can  write 


where 


E a (I)  J_(nir)  . 
m»0  ^ 


We  note  that  the  sequence  {h^  } can  be  treated  as  a matrix  multiplied  by  a 
vector  (in  an  infinite  dimensional  vector  space).  That  is. 


where  a is  a row  vector  of  polynomial  coefficients 

a = tSo*  a, a^j] 

and  B Is  the  (N+l)  x ® transformation  matrix  whose  (m,n)*-^  element  Is  given 
by 

b-n  ■ • 

mn  m 

Observe  that  the  transformation  matrix  B is  independent  of  the  function 
H(x),  and  therefore  a different  transformation  matrix  need  not  be  computed 
for  each  different  H(x).  A procedure  for  the  computation  of  the  terms 
jm(nn)  is  contained  in  section  IV. 

Now  we  consider  the  case  where  H(x)  is  not  a polynomial.  Assume  that 
H(x)  belongs  to  L2[*"ir,ir],  and  therefore  possesses  a Legendre  polynomial 
series  expansion  convergent  In  L2["iT,it].  The  coefficient  vector  becomes 
infinite  dimensional,  and  the  transformation  matrix  B is  doubly  infinite. 

As  a result,  Eq.  (3b)  is  written  as 

h - r a-d)"*  (**) 

It  can  be  shown  [I]  that  the  sum  in  Eq.  (A)  converges  uniformly  in  n.  That 
is,  for  any  e > 0,  there  exists  a K such  that  for  all  M ^ K 

hn  - b„(M)|  < e 

for  all  n,  where 

h (H)  • 1 a (o'"  j_(nir). 

" m-0  ^ ^ . 

Thus  the  sequence  (h-(M)}  is  a uniform  approximation  to  the  sequence  of 
Fourier  series  coefficients  for  the  function  H(x). 

Now,  we  consider  the  inverse  transformation  denoted  by  B~^.  The  com- 
putation of  B'i  directly  from  the  matrix  B can  be  performed;  however,  it  is 
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best  to  compute  the  elements  of  B”*  in  a fashion  similar  to  the  computation 
of  the  elements  of  B.  First,  we  note  that  the  coefficients  {ap}  in  Eq.  (1) 
are  given  by  the  relation  [1] 


a 

n 


2n+l 

2it 


/ HU) 


dx 
n IT 


(6) 


Using  the  Fourier  series  expansion  for  H(x)  found  in  Eq.  (3a),  we  find 


a - (2n+l)  Z h 
n 

IDB-os 


1 r D J 

m 27/ 

-7T  ^ 


Arguing  as  in  the  derivation  of  Eq.  (2),  we  have 


(7) 


If  we  define 


then 


a *=  (2n+l)  Z h [(-i)*’  j„(tmi)]  . 
m»““> 


^mn  “ (2n+l)  (-1)"  j^(nm)  , 


- b"’  ; 


(8) 


(9a) 


(9b) 


that  Is,  the  matrix  elements  of  B“'  are  given  by  b^.  Consequently,  Eq.  (8) 
may  be  rewritten  as 


a ••  h B ' , 


(10) 


where  h is  the  row  vector  of  Fourier  series  coefficients  and  a is  the  row 
vector  of  polynomial  coefficients.  In  practice,  a finite  number  of  elements 
of  h are  computed  by  use  of  the  FFT  algorithm;  we  then  perform  the  vector 
multiplication  indicated  by  Eq.  (10)  to  compute  a finite  number  of  elements 
of  the  vector  a.  This  procedure  for  computing  the  Legendre  polynomial 
coefficients  for  a function  can  provide  a great  reduction  in  computation 
time  as  compared  to  the  direct  evaluation  of  the  integral  in  Eq.  (6). 

In  the  next  section,  we  discuss  a technique  for  the  implementation  of 
nonlinear  operations  by  use  of  a simple  recursive  digital  structure. 

III.  DIGITAL  ZNL  IMPLEMENTATION 


We  begin  with  the  following  recursion  relation  for  Legendre  Polynomials 


(n+1)  (x)  ” (2n+l)  X Pj,(x)  - n P^_|(x)  , (11) 

where 

Pq(x)  - 1 

pj  (x)  “ X , "I  .1  X £ I . 

It  is  straight  forward  to  implement  this  relation  with  a recursive  digital 
network.  Then,  given  this  method  by  which  to  generate  {Pf,(x)  } fur  a given 
X,  and  given  the  coefficients  (a^  )»  the  nonlinear  function  H(x)  is  easily 
constructed  as 

N 

H(irx)  5 E a P (x),  -I  < x < I.  (12) 

n-0  " " " " 

Figure  1 contains  one  possible  implementation  of  Eqs.  (11)  end  (12),  In 
Fig.  1 a recursive  network  is  used  to  generate  Pn(x)  for  successive  values 
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of  n.  With  the  coefficients  {a^}  stored  in  memory,  the  products  and  sum 
indicated  by  Eq.  (12)  are  then  performed.  We  note  that  the  function  H(x) 
impiemented  by  this  network  may  be  easily  changed  simply  by  changing  the 
coefficients  (an);  consequently,  except  for  a memory  unit,  the  network 
remains  unchanged  regardless  of  the  function  H(x). 

Suppose  that  in  order  to  obtain  acceptable  accuracy  at  the  system  out- 
put it  is  necessary,  to  make  N = 50  in  Eq.  (12).  Thus,  given  a value  x we 
must  wait  for  50  terms  to  be  accumulated  before  we  have  the  output  H(x). 
Therefore  this  digital  2NL  could  reduce  the  system  data  rate  by  a factor  of 
50;  however,  for  a price,  it  is  possible  to  eliminate  this  rate  reduction. 
One  approach  would  be  to  simply  use  a faster  clock  rate  with  the  ZNL  . 

structure.  Another  approach  is  to  construct  N = 50  such  2HL  structures;  the 
system' s fi  rst  i nput  sample,  call  it  x^ , goes  into  the  first  ZNL;  the  second 
sample,  X2»  goes  into  the  second  ZNL  and  so  on.  Now,  just  as  x^q  is  going 
into  the  50f*^  ZNL,  we  are  getting  the  output  H(xj)  from  the  first  ZNL. 
Consequently,  after  an  initial  delay  of  50  samples,  50  separate  parallel 
networks  can  operate  at  a rate  50  times  greater  than  a single  network. 
Clearly,  there  are  compromises  that  can  be  made  between  hardware  cost  and 
data  rate. 

Potential  applications  of  this  ZNL  structure  include  adaptive  nonlinear 
systems.  In  such  an  application  the  coefficients  (a^  } can  be  continually 
updated  by  use  of  the  fast  algorithms  previously  discussed. 

IV.  COMPUTATIONAL  CONSIDERATIONS 

In  this  section  we  discuss  two  different  methods  for  the  computation 
of  Jn(n''')«  We  begin  by  noting  the  following  relation  for  generating  values 
of  J^(z)  [2] 

Jn(z)  “ sin  z + (-1)"''’'  f.„_|  (z)  cos  z, 

where 
and 

- (2n*l)  z"’  f„(z)  . 

Consequently, 

j (mw)  - (-l)"'*’*'*^  f ,((mr),  n - 1,2,... 
n •n*i 

where 

il , m ■ 0 
0,  otherwise  . 

The  recursion  in  Eq.  (13)  appears  to  provide  a relatively  simple  method  for 
computing  jp(nm);  unfortunatly , this  relation  is  numerically  unstable.  We 
have  found  that  even  with  double  precision  arithmetic  error  accumulates 
rapidly.  When  using  this  recursion  relation,  we  have  produced  overflow 
conditions  for  n ■ 20  and  z ■>  n,  2ir  on  aUNIVAC  1108  digital  computer;  this 
Is,  of  course,  unacceptable. 

Fortunately,  we  have  developed  a numerical  technique  that  appears  to  be 
very  stable.  First,  we  note  that  the  recursion  relation  for  generating 
Pp(x)  in  Eq.  (II)  is  very  tolerant  of  roundoff  accumulation.  Then,  we 
observe  that  Eq.  (2)  implies  the  following  approximate  relationship. 


(13) 


to 


V ;;  mitt  Siciin 

j (g(  Ml  twtiii  □ 

I MAItOlHICH  □ 

i lOTIfiailOl 


, : . ilM/UglUIIUTY  WtB 
's '^ToTtf  ECWr 
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tn=-j  + 1 


_i 


which  Is  the  discrete  Fourier  transform  relation. 
Inverse  relationship 


N N 


Consequently,  we  have  the 


(I)"  j^(tmrhi 


q=- J + 1 


-1+1  1 

0 <?•••»  0 


Except  for  aliasing  error,  Eq.  (15)  is  accurate.  In  order  to  reduce  alias- 
ing error  It  Is  necessary  to  increase  N;  generally  speaking,  N must  be 
chosen  separately  for  each  case,  depending  on  the  degree  of  the  polynomial 
Pn(*).  Our  procedure  for  computing  jn(nnr)  Is  as  follows:  using  Eq.  (11), 

values  of  P„(^)  are  generated;  then,  using  Eq.  (15),  values  of  jn(mn)are 

computed.  The  values  of  jn(mit)  used  in  the  following  section's  example  have 
been  computed  by  this  method.  We  now  discuss  factors  considered  in  the 
computation  of  the  coefficients  (a^}. 

If  we  assume  that  H(x)  Is  real-valued,  then  Eq.  (8),  used  to  evaluate 
a^,  can  be  simplified  by  noting  that  the  Fourier  coefficients  {h^}  satisfy 
the  relation  h_^  = h*.  Also,  jn(mTr)  satisfies  the  relat  ion  [2] 

Jj^(-itm)  » (-1)"  j^(nm)  . (16) 

Consequently,  Eq.  (8)  may  be  rewritten  as 

7 M/2 

a ■ 2(2n+I)(-l)^  Z j_(mn)  Re  {h_  ) , n-even  and  / 0 
n , n m 

m-I 


~ N/2 
a„  - 2(2n+l)(-l)  ^ E 

nr-1 


j^(mir)  Im  , n-odd. 


where  Re  {h„,}  and  Im  {h„,}  denote  the  real  and  imaginary  parts  of  h^, 
respectively;  and  the  maximum  value  for  f)  is  N,  which  is  the  number  of 
coefficients  {hp,}  that  are  computed.  By  use  of  Eq.  (17) • the  polynomial 
coefficients  {a„ } are  computed  from  the  Fourier  coefficients  {hf,). 

The  following  section  contains  a numerical  example  computed  with  the 
aid  of  the  results  summarized  In  Eqs.  (15)  and  (I7)> 

V.  AN  EXAMPLE 

In  this  section  we  will  demonstrate  the  computation  of  the  polynomial 

coefficients  {a  } for 
n 

H(1TX)  « pj(x)  + P2(x) 

- i (5x^  + 3x^  - 3x-l),  -I  < X < i . 

By  comparison  with  Eq.  (I),  we  see  that  ail  a„  > 0 with  the  exception  of 
*2  " *3  ” example,  ail  computations  are  single  precision 

arithmetic.  The  Fourier  coefficients  th^}  are  computed  by  use  of  the  OFT 


relation: 
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h 


m 


I 

?r 


N/2 

Z 


H&) 


tn 


N 

2 ’ 


and  the  terms  } are  computed  by  use  of  Eq.  (15).  Finally,  the 

coefficients  {ap}  are  computed  by  application  of  Eq.  (17).  These  computa- 
tions are  performed  for  three  different  values  of  N = N:  N = N = 8192, 

1*096,  and  201*8.  Table  I contains  typical  values  of  {a  }. 


“FT 

20^*^[ 

509? 

5T92  ~ 

True  Value 

a 

-4.880  X 10"^ 

-2.441  X 10“^ 

-L 

-1.220  X 10  ^ 

0 

o 

®1 

-1.463  X lo'^ 

-7.320  X lO'** 

-3.661  X lo'** 

0 

«2 

0.9976 

0.9988 

0.9994 

1 

*3 

0.9966 

0.9983 

0.9991 

I 

*4 

-4.385  X lO”^ 

-2.195  X lo'^ 

-1.098  X 10"^ 

0 

•s 

-5.353  X 10‘^ 

-2.681  X lO"^ 

-1.342  X lO"^ 

0 

A 

Table  I.  Selected  Values  of  {a  } Computer  for  N > N. 

n 


As  the  value  of  N Is  increased  the  effect  of  aliasing  is  reduced  in  our 
computation  of  {jn(rnir)  } and  {h^,}.  Aiso,  as  the  vaiue  of  N is  increased,  the 
summation  terms  of  Eq.  (17)  become  more  accurate.  Inspection  of  Table  1 
reveals  that  as  the  value  of  N = 0 is  doubl^ed,  the  error  in  computation  is 
halved.  We  ask,  to  which  parameter,  N or  N,  is  the  error  most  sensitive? 
Some  insight  may  be  obtained  by  examining  the  value  of  a^  in  Table  1;  the 
value  of  aQ  is  independent  of  the  value  of  N (see  Eq.  17).  Because  the 
computation  error  for  ag  is  halved  as  N doubles,  the  implication  is  that 
much  of  the  error  is  due  to  aliasing  and  is  dependent  on  N.  This  |s 
Illustrated  in  Table  2,  where  the  parameter  N is  held  constant  at  N ~ 2048. 


N 

204b 

1*095 

8192 

True  Value 

a 

-4.880  X 10  ^ 

-2.441  X lo"** 

-1.220  X 10"^ 

0 

0 

•1 

-1.463  X 10"^ 

-8.322  X lO"^ 

-6.545  X 10’^ 

0 

*2 

0.9976 

0.9988 

0.9994 

1 

a. 

0.9966 

0.9981 

0.9985 

1 

3 

•u 

-4.385  X 10"^ 

-2.192  X 10"^ 

-1.096  X lO"^ 

0 

•5 

-5.353  X 10"^ 

-3.049  X 10‘^ 

-2.399  X 10"^ 

0 

Table  2.  Selected  Values  of  {a^}  for  N “ 2048. 


First  note  that  In  Tables  1 and  2,  the  values  of  (an)  for  N “ 2048  are 
Identical  because  the  same  value  of  A ■ 2048  is  used^in  both  cases.  Upon 
comparing  the  two  tables,  we  note  that  the  value  of  N does  affect  the  com- 
putation error;  for  some  value  {a^)  the  effect  is  negli gable  while  with 
others  a factor  of  2 difference  is  noted. 

The  results  of  this  section  indicate  the  importance  of  computing  the 
elements  of  the  matrix  very  precisely,  specifically  {jp(mir) }.  If  this 
computation  is  accomplished  with  the  DFT  procedureouti  ined  previously,  then 
a large  value  of  N should  be  chosen  for  the  OFT;  we  emphasize  that  this 
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computation  need  be  made  only  once,  and  the  results  stored  In  memory. 

V.  DISCUSSION  AND  SUMMARY 

In  this  paper,  we  have  developed  a fast  algorithm  for  the  computation 
of  the  Legendre  polynomial  coefficients  of  a function.  Combining  this  fast 
algorl thm  wi th  a digital  implementation  of  the  Legendre  polynomial  recur- 
sion relation,  it  is  possible  to  implement  an  adaptive  nonlinearity.  Gen- 
erally, the  more  complicated  the  non-linearity  to  be  implemented,  the  higher 
will  be  the  degree  of  the  polynosnial  approximation  required;  this  in  turn, 
reduces  the  throughput  rate  of  the  non-linear  device  unless  we  are  willing 
to  pay  a greater  price  in  hardware.  Nevertheless,  by  using  this  approach 
one  acquires  a large  degree  of  flexibility  in  specifying  the  characteristics 
of  the  nonlinearity. 

In  the  past,  most  non-linear  devices  used  in  practice  have  had  very 
simple  characteristics;  in  part,  this  has  been  due  to  the  relative 
difficulty  in  implementing  more  complicated  devices.  However,  as  signal 
processing  schemes  become  increasingly  sophisticated,  the  need  for  complex 
or  adaptive  non-linear  devices  is  growing.  One  very  interesting  application 
Is  In  robust  or  distribution  free  signal  detection  schemes  [3]  where  an 
adaptive,  easily  impimentable,  non-linearity  can  be  of  great  benefit. 

In  addition  to  its  application  to  non-linear  systems,  we  should  not 
overlook  the  advantage  of  having  a fast  algorithm  for  the  computation  of 
the  polynomial  expansion  for  almost  any  function. 

It  should  also  be  noted  that  there  are  many  special  function  pairs 
with  similar  Fourier  transform  relationships  that  may  be  exploited  in  a 
fashion  similar  to  that  presented  in  this  paper  [!»].  The  fact  that  there 
may  not  be  a simple  expression  for  computing  B“'  should  not  be  a deterent 
from  performing  the  computations  because  B"'  need  only  be  computed  once 
and  then  stored  in  memory  for  all  future  requirements. 
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